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ABSTRACT 

We use kinetic particle-in-cell and magnetohydrodynamic simulations supported by 
an observational dataset to investigate magnetic reconnection in clusters of null points 
in space plasma. The magnetic configuration under investigation is driven by fast adi¬ 
abatic flux rope compression that dissipates almost half of the initial magnetic field 
energy. In this phase powerful currents are excited producing secondary instabilities, 
and the system is brought into a state of ‘intermittent turbulence’ within a few ion 
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gyro-periods. Reconnection events are distributed all over the simulation domain and 
energy dissipation is rather volume-filling. Numerous spiral null points interconnected 
via their spines form null lines embedded into magnetic flux ropes; null point pairs 
demonstrate the signatures of torsional spine reconnection. However, energy dissipa¬ 
tion mainly happens in the shear layers formed by adjacent flux ropes with oppositely 
directed currents. In these regions radial null pairs are spontaneously emerging and 
vanishing, associated with electron streams and small-scale current sheets. The number 
of spiral nulls in the simulation outweighs the number of radial nulls by a factor of 
5-10, in accordance with Cluster observations in the Earth’s magnetosheath. Twisted 
magnetic fields with embedded spiral null points might indicate the regions of major 
energy dissipation for future space missions such as Magnetospheric Multiscale Mission 
(MMS). 

Subject headings: magnetic reconnection: null points, simulations: particle-in-cell, 
MHD 


Introduction 


Magnetic reconnection is a fundamental process believed to be the main mechanism of fast 


energy release in magnetized astrophysical plasma (Priest Sz Forbes 2000). It is ignited at kinetic 


scales where ideal conditions for plasma are broken. At the same scales the energy cascade produced 
by turbulence becomes dissipative, hence the two processes: turbulence and magnetic reconnection 
should be mutually connected. 

Many scenarios of magnetic reconnection are attributed to the regions where magnetic field 


vanishes, the null points. Those are found on the Sun (Longcope 2005) where they are believed 


to be the sources of energy release of solar flares (Sweet 1958). In the Earth’s magnetosphere 


null points are accompanied by various instabilities and wave activity patterns (Xiao et al. 2006 
Wendel Sz Adrian||2013 ; Deng et al.pOOO ). Detection of null points requires knowledge of the vector 
magnetic fields implying either polarimetric observations or multi-spacecraft measurements. Such 
measurements are not available in the major part of the solar wind yet. Instead, the tubes of twisted 
magnetic fields, or magnetic flux ropes are found in solar wind (Moldwin et al.|[2000 Janvier et al 


2014). The connection of null points and flux ropes, their relation to magnetic reconnection has 


been unclear, which was a main motive of the present paper. 

Magnetic reconnection in null points has been intensively studied in the framework of magne¬ 
tohydrodynamics (MHD). Classification of the null points is based on the linearization of magnetic 
field: the eigenvalues of the magnetic field gradient define if a null is radial or spiral in three dimen¬ 
sions. Radial and spiral nulls degenerate into X and O points in two dimensions, correspondingly 
(]Lau &: Finn 1990 Parnell et al. 1996). MHD theory of magnetic reconnection in null points has 


been derived by Priest Sz Titov (1996), and reconnection regimes were classified by Priest Sz Pontin 
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(2009). Three-dimensional modeling of magnetic reconnection in null points has been performed 


with MHD (Galsgaard Sz Pontin||20lT| and kinetic (Baumann &: Nordlund[2012 ) codes. Most works 


investigate isolated null points or null pairs, however in situ observations suggest that null points 


tend to concentrate in clusters ( 

Deng et al. 

2009 

Wendel & Adrian 

2013). A study of multiple 

null points in 3D was presented by 

Calsgaard &: Nordlund 

(1997 

), however the addressed magnetic 


configuration was force-free and required artificial external driver to force reconnection, while we 


are interested in spontaneous reconnection events. In the recent MHD analysis of Wyper Sz Pontin 
( 2014b|a ) an initial configuration with a single null point was driven to become unstable. Plasma 
then evolved to a turbulent-like state with null point clusters and magnetic flux ropes governing 
the energy dissipation. 

Recently we have proposed a non-conventional non-equilibrium magnetic field configuration 
where spontaneous magnetic reconnection in null points and flux ropes caused extremely efficient 


energy dissipation (Olshevsky et al. 2013). In the subsequent paper (Olshevsky et al. 2015) we 


have found that plasma was confined in the flux ropes and formed structures similar to Z-pinches. 
Efficient energy dissipation in the system was attributed to the secondary magnetic reconnection 
driven by the instabilities in these pinches. Here we zoom in on the null points in the same con¬ 
figuration, analyze their topology and associated energy dissipation. By comparison of kinetic and 
MHD simulations we have found that the initial relaxation of the system was driven by the large- 
scale fluid modes. These modes released a substantial amount of magnetic energy, but brought 
plasma in the kinetic simulation into a turbulent state within only a few ion gyro-periods. Investi¬ 
gation of different energy dissipation indicators and power spectra allowed drawing out a possible 
mechanism of efficient magnetic energy dissipation in space plasmas. Our idea is that bending and 
interaction of the flux ropes (with embedded currents and spiral nulls) drive secondary instabilities, 
for instance at the shear layers between the adjacent ropes. In such configuration magnetic energy 
cascades from the large (flux rope) to the smallest (electron) scales where electrostatic turbulent 
fluctuations are excited that convert magnetic energy to particle heating. 


2. Simulations 


In this section the results of two simulations are discussed: the main subject of this paper, a 
kinetic PIC simulation described in Appendix [A| and a resistive MHD simulation (Appendix [^. 
Two simulations were carried out on the Cartesian grid with the same number of cells, 400^, and 
with the same initial magnetic/gas energy ratio. The magnetic configuration under study is fully 
periodic, magnetic field is given by 
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The simulation boxes had dimensions Lx = Ly = Lz = 20 in the corresponding code units (ion 
inertial length di in the PIC simulation, Appendix [A]) . Although the adimensional units in the 
MHD code may be scaled to any physical range, it is inappropriate to scale them to di, because 
the essence of the MHD model is its large (comparing to, e.g., ion) scales (Appendix [B|) . Therefore 
we don’t directly compare times and distances between the two runs. However, some relative 
quantities such as the kinetic/magnetic energy ratio can be quantitatively compared. We have set 
the initial magnetic field amplitude to Hq = 1 in the MHD simulation, and to Bq = 0.0127 in the 
PIC simulation to match the initial kinetic/magnetic energy ratio. Initial uniform density was set 
to unity in each code’s units. 

The magnetic field configuration described by Equations[^is divergence-free, but not force-free. 
Initial density distribution is uniform, hence forces at t = 0 are not balanced. Pressure imbalance 
triggers explosive relaxation during which about a half of the initial magnetic energy is released 
(Figure]^ a, b). In addition to energy, magnetic field and density distributions are also comparable 
during this phase as seen from comparison of panels c and e of Figure Both observations suggest 
that the initial evolution is driven by large-scale MHD modes. 


We use the Poincare index method, as described in Appendix [Cl to detect and classify magnetic 
null points in the simulations. The classification that we use ( [Cowley 1973; Lau Sz Finn 1990) is 
based on the VB eigenvalues, which should sum up to zero due to the divergence-free condition. 
When all three eigenvalues are real the null is called radial. Radial nulls degenerate into X points 
in 2D. When two eigenvalues are complex, the null is called spiral. The topology of magnetic field 
in the vicinity of a spiral null resembles a plasmoid or a helical flux rope; a spiral null degenerates 
into an O point in 2D. Further division of nulls is defined by the signatures of the real parts of 
the VB eigenvalues: when two eigenvalues are positive, the null is of B type, also called positive 
(Longcope 2005); when two eigenvalues are negative, the null is of A type {negative)] similarly, 
spiral nulls are divided into As and Bs types. Our initial magnetic configuration contains lines 
composed of essentially two-dimensional, non-generic O points that are structurally unstable in 3D 
( Greene[[1988 ). The O points evolve into As and Bs nulls just after the beginning of the simulation, 
and the 3D null lines consist of interconnected spiral nulls As and Bs (Fig. [^c, e). 


The energetics of the two simulations is illustrated in panels (a, b) of Figure The time scale 
in the MHD run is arbitrary, and the reader should not be confused by the similarity of energy 
plots: only the initial phases are similar. Panels (c-f) of Figure display the snapshots of the 
PIC (c, d) and MHD (e, f) simulations taken at two moments: during the explosive relaxation 
and at a later stage (marked by the vertical lines in panels (a, b). In the beginning the null lines 
are fenced in the helical magnetic field lines. Because the initial density is uniform, magnetic 
tension squeezes plasma towards the null lines where high-density regions are formed, as illustrated 
by the grey shade in panels (c, e) of Figure (only 3 out of 9 null lines are shown for clarity). 
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The process of compression is followed by the reverse expansion, resulting in pulsations visible in 
kinetic and magnetic energy plots in panels (a, b) during some one-third of the total duration of 
the simulations. Notably, the major part of magnetic energy released during the first compression 
in the PIC simulation is transferred to ions. These ions create powerful currents along the null 
lines, and their collective behavior is well approximated by the fluid model. 


When the pressure imbalance is compensated, the behavior of plasma diversifies in the MHD 
and PIC simulations. Such dissimilarity is caused by the intrinsic discrepancy between the spatio- 
temporal scales considered by the two physical models. In the kinetic simulation there is no 
numerical resistivity, and the currents along the null lines do not dissipate for relatively long 
period. Hence the plasma confinement holds, and its behavior is dominated by the current channels 
( ]Qlshevsky et al.]|2015 ). In the MHD simulation the currents are ruined by the resistive dissipation 
and become disrupted quickly. In the absence of a large-scale driver, the simulation becomes 
very chaotic. Figure illustrates the flows in the simulation domain at the same moment as 
in Fig. [^. They indeed show no large-scale coherent structures, with vortices of different sizes 
formed everywhere. When the pressure imbalance is compensated in the MHD run, the spectrum 
of magnetic fluctuations (at well resolved scales) gets a power-law shape with exponent « —1.7 
(Fig- [^), which corresponds quite well to the measured fluid-scale spectrum in the solar wind 


(Alexandrova et al.||2012 


Safrankova et al. 


2013). The fluctuations at scales of a few grid cells are 


damped by the numerical scheme (slope limiter, see Appendix]^, and the spectrum becomes much 
steeper. 


The situation is quite different in the PIC simulation, where specific scales are defined by 
ions and electrons. The absence of numerical resistivity allows currents to hold for relatively 
long periods, and the confinement holds for tens of ion gyro-periods. As shown by [Qlshevsk^ 


et al. (2015), Z-pinches are formed along the sinusoidal null lines along which the currents are 


streaming through twisted magnetic fields. These pinches play major role in the energetics of the 
PIC simulation. Their configuration resembles plasmoids extended to three dimensions ( Markidis| 
et al.||2012 Vapirev et al.||2013) and flux ropes observed in the solar wind (see Fig. and Fig. [^. 


Adiabatic compression of the flux ropes driven by the pressure imbalance releases almost half of 
the initial magnetic energy, but brings the system to a state with a well-defined energy cascade 
(Section]^ and numerous small-scale reconnection events governing energy dissipation (Section]^, 
within only a few ion gyro-periods. Null points are important actors in these processes, and it 
is appropriate to first zoom in, and understand what is happening in nulls during reconnection. 
Since magnetic reconnection is a kinetic process, we will only consider the PIC simulation in the 
remainder of the paper. 


3. Observations 

Of major importance is to understand the relative roles of spiral and radial null points in the 
processes of energy dissipation in space and solar plasmas. Certain indications that spiral null 
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points might dominate over the radial ones were found in observations (Xiao et al. 2006; Wendel 


&: Adrian 2013 Deng et al. 2009) and simulations (Wyper &: Pontin 2014b). Unfortunately, no 


systematic observational survey has been performed yet: main limitation of the widely employed 
Poincare index method of null detection (Appendix is its inability to detect nulls outside the 


spacecraft tetrahedron. The more advanced method proposed by Fu et al. (2015) will provide means 


to perform such a survey (Eriksson et al. 2015), which we briefly demonstrate in this section 


To fill in the gap and put our study into context, we have applied the two methods, Poincare 
index and Taylor expansion, to a set of full resolution (~ 66 Hz) DC magnetic field data from 


the Flux Gate Magnetometer (FGM) (Balogh et al. 2001) aboard Cluster spacecraft (Escoubet 


et al. 2001). The data was taken on 27 March 2002 from 09:45 to 11 UT, when the four Cluster 


spacecraft were in the magnetosheath, downstream of a quasi-parallel bow shock. The same period 


was studied before analysed, in particular, by Retino et al. (2007); Wendel &: Adrian (2013) 


Figure]^ shows an example of null point detection performed over a 3 minute set of magnetic 
field data (Panels a-c) on 27 March 2002, around 10:06 UT — 10:09 UT. Positions of the nulls 
are plotted in panels d-f (X, Y, Z in GSM coordinates, respectively). Null points are located with 
a good accuracy in weak |i?| regions, as indicated by error constraint conditions marked by black 
lines (Appendix [C|) . Notably, positive and negative nulls are commonly found in pairs as visible in 
Figure]^ 10:06:05 UT, 10:07:58 UT, which is also common in our PIC simulations (see below). 

A zoomed view around 10:08:05 UT (Figure]^ g-1) provides an example of transformation of a 
radial A null into a spiral As null and back. At the present stage it is difficult to interpret this result 
as a clear observation of a single null point conversion. The detection procedure assumes linearity 
of magnetic in the vicinity of a null point, whereas magnetic reconnection simulations reveal the 


multiscale picture of reconnection with electron beams, plasmoids, thin current layers (Karimabadi 


et al. 2014), which can significantly affect remote detection of null points in satellite data. 


The numbers of null points of different types are summarized in Table The Taylor expansion 
method has found 443 time steps with nulls, while the Poincare index has found only 64, illustrating 
the aforementioned limitation of the second method. Both methods have detected the dominating 
(80-90%) number of spiral nulls. As and Bs. There was a slight tendency to have more negative 
(A and As) nulls than positive (B and Bs) nulls. 

We used Poincare index method as described in Appendix [C| to find null points in the simu¬ 
lation domain, the results are shown in Figure Our initial magnetic field configuration contains 
8 null points of radial A and B types and (mathematically) infinite number of degenerate 0-type 
nulls. The 0-points, initially forming the null lines, evolve into spiral As and Bs nulls as the sim¬ 
ulation begins, which explains the sharp drop of the ratio of spiral/radial null numbers. In such 
a fashion we artificially create a relaxing system where the spiral nulls outweigh the radial nulls, 
and the spiral/radial null count ratio keeps high (5-10) throughout the simulation, in accordance 


with observations. Wyper & Pontin (2014a) concluded the importance of the spiral null points and 


magnetic flux ropes for magnetic reconnection in their 3D MHD simulations initiated in a com- 
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pletely different magnetic field configuration with a single isolated radial null; clusters of secondary 
null points were forming during the evolution of such a system. Therefore, in addition to the ob¬ 
servational survey, a numerical investigation should be performed in different field configurations 
to either support or drop this conclusion. 


4. Spiral nulls 


Previous works have indicated that energy conversion in the discussed configuration is asso¬ 
ciated with spiral null points connected into null lines embedded into twisted magnetic flux ropes 
( jOlshevsky et ah 2013, 2015[ ). Figure shows magnetic topology and electron streamlines associ¬ 
ated with two adjacent flux ropes with null lines. Positive and negative (Bs and As) nulls alternate 
along the null line. The black magnetic field lines in Figure belong to the fan surfaces of the 
corresponding nulls. The separators of the adjacent nulls are created by intersecting fans; helical 
field lines also connect the nulls with the external magnetic field. This picture is very similar to 


the one used by Wyper &: Pontin (2014a) to describe the ‘secondary bifurcations’ that lead to the 


creation of spiral null pairs in their MHD simulations. Similar topology was observed in a pair of 


nulls detected by Cluster in the Earth’s magnetosheath ( 

Wendel & Adrian 

2013). Reconstruction 

of magnetic topology of another spiral null pair by 

Deng et al. 

(2009 

) revealed a large angle between 


the spines of the two nulls, and a fan-fan separator line. 

Although we haven’t performed a dedicated topological analysis of our simulations, the helical 
field lines surrounding the nulls in Figure a,b suggest their spines to be tangential to the current 
wires. MHD theory would classify magnetic reconnection in such configuration as torsional spine 


reconnection: the currents accumulate along the spines and are co-aligned with them (Priest & 


Titov||19^ Priest &: Pontin|2009 ). In panels (a, b) of Figure]^ a zoom-in on the magnetic topology 
of the null lines reveals a ‘knot’ of field lines in the central flux rope (highlighted by a blue box in 
panel b). No null points are detected in this region, which might indicate an emerging or decaying 
As-Bs null pair. Interestingly, not much energy conversion is associated with this region or the 
majority of null points directly. The divergence of the Poynting vector V • S and the work of the 
field E-J on particles are high in the regions where adjacent current channels bend or interact with 
each other, and where clusters of nulls are found. This observation is supported by other indicators 
such as electron heating or the violation of frozen-in condition already presented in Figure 4 by 


Olshevsky et al. (2015). Note, that although the maximum absolute values of V • S is substantially 
higher than the maximum of E • J, domain-averaged value of the latter quantity is about 10 times 
larger than the average value of Poynting vector. 


The in-plane electron velocities v^x and Vez in the slice through the simulation domain show 
several shear layers where the flux ropes interact or bend (Figure]^ c, d). Shear layer indicated by 
the blue and green electron velocity vectors (same as in panels a, b), and the symmetrical one in 
the upper right corner of the slice demonstrate oscillatory pattern in v^z- The wavelength of the 
oscillation is of the order of di, however, our preliminary simulations with higher resolution show 


























that these oscillations have fine structure, and they will be addressed in detail in a subsequent 
paper. 


5. Radial nulls 


Figurej^presents a typical A-B null pair spontaneously emerged in the simulation at Qcit = 42.5 
on the interface of the two flux ropes (red and yellow field lines in panel a). The distance between 
the nulls is Id*, they disappear in the subsequent simulation snapshot, hence their lifetime is no 
longer than 712^^. The interconnection between these radial nulls is more complex than in the 
previously shown case of spiral nulls. The field is non-linear, field lines forming the fans and spines 
of the nulls are bent and twisted, and reconstruction of such topology from observations would 
require more sophisticated methods. Magnetic held lines enter the A null (cyan) in the curved fan 
surface, and exit it along the spine line, and visa versa for the B null (blue). The fan of the A 
null and the spine of the B null are formed by the same group of held lines that start from the X 
domain boundary (in the bottom of panel a). On the bottom, these held lines encircle the yellow 
hux rope, while upper they twist around the red one, joining the topologies of the two ropes. The 
short lifetime of the nulls, their close location, simultaneous appearance and vanishing suggest they 
are connected by a separator, which, most likely, lies at the intersection of the fan surfaces of the 
two nulls. 


A stream of electrons co-aligned with the spine of the positive (B) null (black arrows in Fig¬ 
ure §b) approaches the null point and rehects towards its negative (A) companion. After passing 
the A-point electrons scatter and form a small-scale current sheet surrounded by the twisted mag¬ 
netic held lines. This process is supported by the positive work of the held on the particles and 
generation of EM huctuations in the region where the divergence of Poynting vector is negative. 
The described null pair is the hrst among several A-B pairs created at a later time (see Figure]^. 
Notably, not always null points are born and die in pairs in our simulations and observations: at 
certain moments the number of positive nulls deviate from the number of negative nulls. Theory 
predicts conservation of the topological degree, i.e., each new or disappeared positive null must 
have its negative counterpart and vice versa. In practice errors occur on the stage of null detec¬ 
tion due to the noise in magnetic held measurements and hnite resolution of the numerical grid 
or spacecraft instruments. Additional errors are introduced during null classihcation, especially 
when the eigenvalues of the VB have very small imaginary (real) parts. Inhuence of these errors 
on the observations of null points are addressed in Fu et al. ( |2015 ); priksson et al. (2015), see also 
Appendix 
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6. Energy cascade 


As we pointed out in previous sections, energy dissipation in our simulation is associated 
with electron and ion beams streaming through null points, small-scale curent sheets and sheared 
motions. Especially important the large-scale flows are in the first phase of “implosive” relaxation 
when the energy of ion beams is about a half of the magnetic field energy (see Fig. a). Two- 
stream and shearing instabilities associated with these flows are expected to drive electrostatic 
turbulent fluctuations and lead to creation of energy cascade. A somewhat better resolved 3D 


kinetic simulation ( 

Daughton et al. 

2011 

) of flux ropes produced by magnetic reconnection in 

a Harris sheet (Harris 

1962 

) has indeed shown ‘intermittent’ turbulence, analyzed in detail by 

Leonardis et al. ( 

2013) 

. Of course, our simulations do not have external driver which is typically 


present in turbulent studies, and, at our best, we can attribute it as a case of intermitent turbulence. 

Indeed, in a few ion gyro-periods magnetic energy spectrum obtains a characteristic shape 
(red curve in Fig. a) with four distinct breaks. The first break at kdi = 1 is a transition from 


fluid to kinetic (ion) scales, similar to the one observed in solar wind (Alexandrova et al. 2012 ). 
Unfortunately our simulation box is too small to properly represent the fluid range in the spectrum. 
The second break is associated with electron inertial scales kde = 27r and also corresponds to the 
one detected in solar wind by Sahraoui et al. (2013). The magnetic energy spectrum between 


kdi = 1 and kdg = 1 is rather smooth power-law with a best-fit exponent 7 ps — 5 throughout 


the simulation. Such steep shape is found at the sub-electron scales in solar wind (Sahraoui et al 


2013). 


The third break in the B spectrum kdi ~ 50 makes it milder, while the fourth steepens it 
further, in discordance with the measurements of Sahraoui et al. (2013). These scales, however. 


correspond to the wavelengths shorter than 2 grid points, where the implicit scheme is damping 
numerical noise. Panel b of Figure shows the spectra of electric field fluctuations. To damp 
numerical noise and stabilize the solution, smoothing is applied to the electric field derived at each 
time step. The smoothing operates at five-point stencil, and small-scale E fluctuations are damped 
to machine precision. The corresponding break is seen in the electric field spectra at kdi ~ 20. 
Although smoothing is applied only to electric field, which energy budget is substantially smaller 
than that of the magnetic field, it might also affect the shape of the magnetic spectrum at the 
smallest scales. 

The distinct features of the electric energy spectrum at Dcjt = 3.5 are the hump at ion inertial 
scale, followed by a local minimum at kdi = 27r, and another hump at kre = 27r (electron gyroradius 
Te is computed from the domain-averaged average magnetic field and electron thermal speed at 
ti^ci = 70.9). The first hump disappears at later time moments, and we associate it to the initial 
compression of the magnetic flux ropes, excitation of ion currents and acceleration of ions to non- 


thermal speeds (Olshevsky et al. 2015). The peak at electron scales that may be blurred due to 


the aforementioned smoothing. 

The spectra of electron velocity fluctuations (Figure c) have a well defined peak at the 
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smallest, sub-electron scales at all time moments. The energy contained in these scales increases 
with time, indicating heating of electrons due to the dissipation of magnetic energy at the smallest 
scales. Notably, the power of velocity fluctuations and electric field fluctuations start increasing at 
frequency kdi = 10. The higher frequency fluctuations of electric field, however, are damped by 
numerical smoothing. The variation of the velocity spectra with time allows to conclude that in 
the beginning (red) most of the energy is contained in the current channels with diameters of about 
di; later (green and blue), the energy redistributes from the large scales to the small scales, with 
a distinct local minimum between ion and electron scales. In the absence of (externally driven) 
large-scale fluid motions the energy from ion-scale magnetic structures is redistributed to electron 
scales, where electron velocity fluctuations are excited. 


7. Summary 

We have investigated magnetic reconnection in a cluster of null points by means of kinetic 
particle-in-cell simulation. The simulation is initiated from an unbalanced configuration, abrupt 
relaxation of which excites powerful currents and brings the system to a turbulent state within a 
few ion gyro-periods. Major role during relaxation is played by adiabatic compression of magnetic 
flux ropes: it excites powerful currents along the embedded null lines, and accelerates ions to 
suprathermal speeds. Large-scale fluid oscillations are excited at this stage, as confirmed by an 
accompanying MHD simulation performed in the same magnetic configuration. 

When the pressure imbalance is compensated, confinement is ruined in the MHD run, and 
the system becomes very turbulent. Energy spectra obtains a distinctive power-law spectra with 
exponent ~ —1.7 in the well-resolved range, becoming steeper at the shorter wavelengths where the 
dissipation is caused by numerical scheme. Multiple null points are created all over the simulation 
domain, but an MHD simulation doesn’t allow for an insight into the micro physics of magnetic 
reconnection in turbulence. 

In the PIC simulation distinct spatial scales are dictated by particles, plasma confinement 
is preserved, and the dynamics in the stationary dissipation phase is dominated by interacting 
magnetic flux ropes. The flux ropes in our simulation are formed along null lines, the sequences 
of spiral null points connected via their spines. Magnetic reconnection in these null points can be 
classified as spine reconnection, with energy being transferred to particles mainly in the regions 
of flux rope bending or interaction. On the interfaces of oppositely directed adjacent flux ropes 
the shear layers are created that display oscillatory pattern with a wavelength of the order of ion 
inertial length. A simulation with higher spatial resolution is required to analyze this instability. 

The number of null points in the simulation varies with time, with null pairs spontaneously 
emerging and vanishing. A typical pair of radial null points is created on the interface of the two 
interacting flux ropes. Such pair is usually short-lived, it reconnects in a few ion gyro-periods. The 
distance between the nulls in the newborn pair is about one ion inertial length. Visualization of the 
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local magnetic topology suggests that the separator of the pair lies in the fan surfaces of the nulls. 
Magnetic reconnection in such pair is characterized by a stream of electrons entering the system 
along the spine of the positive null, refracting towards the negative null, and then scattering away 
in a small-scale current sheet. 


Magnetic reconnection changes the number of null points in the simulation domain, however 
there always are 5-10 times more spiral than radial nulls. We have applied two methods of null point 
detection: Poincare index and Taylor expansion, to a one-hour period of Cluster measurements in 
the Earth’s magnetosheath. The second method has found seven times more nulls (443) than the 
hrst one (64), and in both cases more than 80% of nulls were spiral. For a spiral null pair from the 
same period of Cluster observations Wendel &: Adrian (2013) have deduced magnetic field topology 
and revealed a spine-aligned current, which corresponds well to our simulations. 


Our kinetic simulations, supported by the results of in situ observations in the magnetosheath 
suggest a new mechanism of magnetic energy dissipation in turbulent space plasma. It is governed 
by bending and interaction of magnetic flux ropes that produce the cascade of magnetic energy 
to the sub-electron scale, where small-scale electron fluctuations are excited. Electron velocity 
fluctuation spectra gets a characteristic maximum at sub-electron scales. Locally, the dissipation 
is attributed with clusters of interconnected spiral null points with co-aligned current ropes, short¬ 
living pairs of radial nulls, shear layers and small-scale current sheets. 


Authors are thankful to Yu. Khotyaintsev, A. Vaivads, and H. Fu for useful discussions. This 
research has received funding from the European Commission’s FP7 Program with the grant agree¬ 
ment eHeroes (project 284461, www.eheroes.eu). E.E. and S.M are supported by the Swedish VR 
grant D621-2013-4309. The simulations were conducted on the computational resources provided 
by the PRACE Tier-0 projects 2011050747 (Curie) and 2013091928 (SuperMUC). 


A. PIC simulation 


The kinetic simulation of collisionless plasma was carried out using semi-implicit fully electro¬ 
magnetic PIC code iPic3D (Markidis et ah 2010). The code solves the time-dependent Maxwell 
equations for the fields given on a stationary grid, and the equations of motion for computational 
particles derived from Vlasov equation. Each computational particle represents a blob of real par¬ 
ticles (ions and electrons) that are close to each other in 6D phase space. The physical units in 
the code are normalized to the corresponding plasma parameters: proton inertial length di, proton 
plasma frequency iVpi, and proton mass m*, hence magnetic field unit is ruiUJpi/e. 


The simulation was carried out in a cubic box of size 20 di and the resolution of 400 cells in each 
dimension. Two species, ions and electrons were considered, with mass ratio mj/rue = 25, and 64 
particle of each specie per cell. The time step was set to 0.15ti;pi, satisfying the finite-grid stability 
criterion; the total duration of the run was 70 ion gyro periods Particles were initiated with 
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a Maxwellian distribution with thermal speed in each dimension uth,e = 0.02 for electrons, and 
uth,i = 0.0089 for ions, which corresponds to the temperature ratio Tj/Tg = 5 typical in the Earth’s 
magnetosheath plasma. 


Initial domain-integrated magnetic/thermal energy ratio was set to Wmag/Wth = 1-38 and 
decreased to 0.07 in the end of the simulation, passing through the stages of different plasma j3. 
Given these parameters, the electron inertial scales were resolved substantially well, however to 
resolve electron gyroradius at all stages and to study the fine structure of induced instabilities, 
higher resolution simulations or mesh refinement techniques (Innocenti et al. 2015) are required. 


B. MHD simulation 


A total variation diminishing Lax-Friedrichs (TVDLF by [Toth &: Odstrcil ( ]1996 )) scheme 
with explicit time stepping (second-order accurate in space and time) was implemented in the code 
infrastructure of iPic3D (hereinafter referred to as ‘iPic3D-MHD’) based on the serial MHD code 
by Kiehas et al. (2006). 

iPic3D-MHD solves conventional resistive MHD equations cast in dimensionless form (see, e.g.. 


Ma & Bhattacharjee (2001)). The ideal MHD equations were used for the present study 


Ft 


de 

Ft 


+ V 


| + v.(,o = 0. 


2 


ve + V [p + 

9B 


/ B^\ 



= 0 

-BBu 1 

= 0 


dt 


+ V • (uB — Bu) = 0, 
2 ^2 


_ p . pv- 

where e is the total energy density, 7 = 5/3 is the specific heats ratio. The energy components shown 
in Figurej^ are given by: magnetic field energy Wb = 5 / B^dV, thermal energy Wth = 7 ^ / pdV, 
bulk flow energy W^uik = 5 / pv'^dV. 


The numerical diffusivity required for the stability of the simulation was provided by TVD 
slope limiters (Toth fc: Odstrcil]|1996). We have tested all three slope limiters proposed by Toth & 


Odstrcil (1996), and haven’t found a notable influence on the energetics of the simulation. 


Magnetic field at time t = 0 was set according to Equations using Hq = 1 in dimensionless 
units. Plasma density p = 1 and pressure p were set uniform at t = 0. Initial plasma pres- 
snre p equals to pi + Pe^ the sum of the initial ion and electron pressures in the PIC simulation 
{{me/mi)uff^ ^ + u^hi)/BQ = 0.6). The integrated initial magnetic to thermal energy ratio was 
Ws/Wth = 1.38, same as in the PIC simulation. 
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The divergence-free condition for magnetic field (VB = 0) is satisfied at t=0. However, the B 
field can slowly accumulate magnetic charges over the course of the run, unless corrected. Numerical 


lioo| ). 

A uniform Cartesian grid (400 x 400 x 400 cells) with equal spacings in all three dimensions 
was used. The number of grid points was identical to that in the PIC simulation (Appendix A) to 
ease the comparison. 


magnetic field divergence is cleaned each 10*^ time step by means of the projection method (Toth 


C. Location and classification of null points 


The problem of locating a magnetic null is essentially a problem of finding a root of a continuous 
divergence-free vector field. The widely adopted method of null point location is based on the 


topological degree, or Poincare index, and was introduced by Greene (1992). A topological degree 


estimates a number of roots (nulls) in a closed region of space; it is non-zero when an odd number 
of roots are enclosed in the region. Commonly used assumption is that the separation between 
spacecraft is so small that the non-zero topological degree corresponds to exactly one null. 


The method introduced by Greene (1992) used a rectangular box as a volume for which the 


topological degree was computed. However, since Cluster mission consists of four spacecraft, the 
method was adjusted to find nulls in a tetrahedron. Null detection in the simulation data was 
also performed with the same method for consistency. Each cubic mesh cell was divided into hve 
tetrahedrons with vertices at the mesh nodes. 

Computation of the topological degree is based on finding a solid angle between vectors in 3D. 
Unlike other implementations we used the formula for solid angle enclosed by 3 vectors proposed by 


van Oosterom & Strackee (1983); it appears faster and more stable than the traditional implemen¬ 


tation based on the Cosine theorem. In particular, there is no need for zero-denominator checks 
when using a present-day programming environment: errors are handled by the arctan2 function. 

The magnetic field is assumed to be linear in the vicinity of a null point. Linear expansion 


of the magnetic held is a basis for the classihcation of nulls introduced by Cowley (1973); Greene 


(1988); Lau & Finn (1990). The eigenvalues of the VB tensor dehne the topological type of the 


null. We detect the null point type by estimating the VB as described in Khurana et al. (1996) 


Obviously, a shortcoming of the Poincare index method is that it can only detect a null point 
inside the tetrahedron. This is hne in the models where the “spacecraft” cover the whole volume 
of the simulation domain, but is a severe limitation in real observations. Therefore, to the latter, 
we have also applied a more advanced method based on Taylor expansion of the magnetic held in 
the vicinity of a null ( Fu et al^|2015 ): 


B (r) = VB (r - tq) , 


(Cl) 
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where r is the location in space, and fq is the location of the null. Inversion of this linear expansion 
gives the null position fq. In general, the equation will always give a position of a magnetic null. 
However, we regard the null reliably identified only if it is located in a rectangular box defined by 
the spacecraft location. The edges of the box in each direction (x, y, z) are given by the maximum 
and minimum positions of all Cluster spacecraft. Only the magnetic null positions found within 
this box are type identified and further evaluated. 


We have used two error constraints to estimate the accuracy of null point detection: 


V B 


Al + A 2 + A 3 

max(VB) 


max{\real{X)\) 


(C2) 


Possible errors of null detection via Taylor expansion are addressed in Fu et al. (2015). For the 
purpose of the current work, we have tested four thresholds set for both error conditions (i.e., when 
at least one error condition is above the threshold, the null is discarded): 10, 20, 30, and 40%. As 
expected, we found no notable difference in the ratio of the spiral/radial null numbers: indeed, the 
threshold only limits the number of the selected data points, not their type. Table presents the 
results for the error threshold of 40%. 
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Fig. 1.— Comparison of PIC and MHD simulations, a, b) Evolution of different components of 
energy in PIC (a) and MHD (b) runs. Domain-integrated energy components are marked: IPs - 
magnetic field; Wth, Wth,e-, Wth,i - thermal energy in the MHD run, and thermal energies of the two 
species in the PIC run; Wbuik, Wbuik,e, Wbuik,i - same for the bulk kinetic energy. Vertical dashed 
lines mark the moments at which the snapshots are made, c, d) Snapshots from PIC simulation 
corresponding to the two moments marked in (a). Null-points are marked with spheres which colors 
designate the null point types. Green magnetic field lines are associated with a null line; magenta 
field lines show the topology of radial nulls. Mass density is shown with shades of grey, e, f) 
Snapshots from MHD simulation corresponding to the two moments marked in (b). Meaning of 
elements is the same as in (c, d). Note, in (c, d, e, f) the regions adgacent to Y boundaries are 
masked for clarity, due to the symmetry of the configuration they don’t carry additional information. 
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k-L 


Fig. 2.— Multi-scale chaotic flows in the MHD simulation, a) Same time moment as in Figure [^f 
shown under a different angle. Same magnetic topology illustration, but the null points made 
smaller for clarity. The T = 10 plane is colored with velocity amplitude, b) Power spectral density 
(PSD) of magnetic field at three time moments in the MHD simulation. The color of the spectra 
indicates the time step index i = 500 (red), i = 5500 (green), and i = 8500 (bine). The vertical line 
indicates the wavenumber corresponding to the wavelength of 5 grid cells. The best-fit exponents 
are qnite far from the expected values in MHD turbulence; the dissipation is purely numerical. 
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Fig. 3.— Detection of null point in Cluster spacecraft data, (a-c) The X, Y, Z components of 
magnetic field measured by FGM. (d-f) The X, Y, Z coordinates of the null points found within 
error margins (black lines), as computed from equation 1. The scmin and scmax label the minimum 
and maximum coordinates of the spacecraft along the corresponding axis (mark the edges of the 
rectangular box of the spacecraft as explained in Appendix [C]) . (g-l) Same as (a-f), but for 0.5 
second interval around 10:08 UT marked by yellow bar in (a-f). 
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Fig. 4.— Evolution of the number of null points of different types in the PIC simulation. The 
colors correspond to the null point color legend in Fig. Two-dimensional 0-type nulls exist only 
in the initial state, often misinterpreted as spiral nulls, the reason for the very high spiral/radial 
ratio at t = 0. 
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Fig. 5.— Magnetic topology of the spiral nulls at Qcit = 35. a, b) A zoom-in on two adjacent 
flux ropes. Spiraling grey magnetic field reveal the fans of the As and Bs nulls connected via their 
spines. Notable is a ‘knot’ in the spiral highlighted by a blue box in (b) that does not contain 
null points. It may be associated with birth or decay of a null pair. Green and blue arrows mark 
electron streamlines in the two oppositely directed current channels. Electrons travel along null 
lines for 10-20 di, and then diverge. The red and blue shade is volume rendering of the divergence of 
Poynting vector (a) and the magnitude of E • J (b). The most intensive energy conversion happens 
in the region where the two flux ropes interact or bend, c, d) A slice through the simulation domain 
at F = 10 di shows the in-plane electron velocity components Vex (c) and Vez (d) with a grey-red 
pallette. Null points, velocity vectors and field lines are the same as in the upper panels. In the 
lower left corner a shear layer is formed between two current channels. On the interface a wave-like 
pattern is visible in (d) that is also seen in the upper-right corner where there is a symmetrical 
shear layer. 
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Fig. 6.— Magnetic topology of a radial null pair at Qcit = 42.5. a) The A-B null pair (shown with 
blue and cyan field lines) is formed on the interface between twisted fields surrounding adjacent 
null lines (marked with red and yellow field lines). The arrows are magnetic field vectors along the 
plotted field lines, b) Same null pair viewed under slightly different angle. Blue and cyan field lines 
depict the topology of the A and B nulls. Black arrows mark electron velocity vectors: electron 
stream enters the B null, bends towards the A null, and scatters away forming a small current 
sheet surrounded by twisted magnetic field lines. Red and blue shade marks the regions where the 
divergence of Poynting vector reaches the highest values. As in Figure these regions correspond 
to the regions of high E • J (with the opposite sign; not shown here). 
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Fig. 7.— Spectra of magnetic field (a), electric field (b) and electron velocity (c) fluctuations at 
three moments in the simulation. The color of the spectra indicates the time moment Qcit = 3.5 
(red), hlcii = 35 (green), and ^Icit = 71 (blue). The vertical lines indicate kde = 27r (dashed), 
kdi = 27r (solid), and kr^. = 2 tt (dotted) computed at the end of the simulation basing on the 
electron thermal speed and domain-averaged magnetic field value. The energy from large-scale 
magnetic features is transferred into small-scale motions and heat (maximum on the right of the 
panel c). To stabilize the numerical solution, electric field is smoothed beyond electron gyroscales 
(sharp drops in panels a, b around the dashed line). 
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Nulls found 

A(%) 

B (%) 

As (%) 

Bs (%) 

Poincare index 

64 

8 

1 

55 

36 

Taylor expansion 

443 

14 

8 

42 

36 


Table 1: Percentage of magnetic nulls of different types detected with Poincare index and Taylor 
Expansion methods in Cluster measurements in Earth’s magnetosheath. 



